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Abstract. The ensemble of random Markov matrices is introduced as a set of Markov 
or stochastic matrices with the maximal Shannon entropy. The statistical properties of 
the stationary distribution tt, the average entropy growth rate h and the second largest 
eigenvalue ly across the ensemble are studied. It is shown and heuristically proven that 
the entropy growth-rate and second largest eigenvalue of Markov matrices scale in 
average with dimension of matrices d as h ^ \og{0{d)) and \v\ ~ respectively, 
yielding the asymptotic relation hr^ '-^ \j2 between entropy h and correlation decay 
time Tc = — 1/log \ v\ . Additionally, the correlation between h and and is analysed 
and is decreasing with increasing dimension d. 
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1. Introduction 

In the information tlieory and matliematical modeling of pliysical processes we often 
stumble upon Markov cliains [1] and Markov or stocliastic matrices |2] , which determine 
the evolution of the first. Let us assume that a Markov chain is based on the set of 
states S = {si}f^i and P{sj\si) is the conditional probability for a transition from state 
Si to Sj then the corresponding Markov matrix M E W^'^ is a collection of conditional 
probabilities 

d 

M,,, = P{s,\s,)>0, J]M,,, = 1, 

where d the dimension of the Markov matrix. Notice that the sum of elements in each 
rows is normalized to 1. The applications of Markov matrices as well as the ways how 
to constuct them are very diverse. In particular interesting is their use in dynamical 
systems for a probabilistic description of the dynamics. For a general introduction in 
this direction see e.g. [3]. For example let us consider a discrete dynamical system 
0* : X — ^ X, where t G N, with the phase space X and the invariant measure /z. By 
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choosing disjoint subsets of phase space {Xi C X : Xj n Xj = for z 7^ j}, which satisfy 
IJ. Xj = X, the Markov matrix M = [Mij]fj^i corresponding to the dynamical system 
can be defined as 

_ /i(0mnx,) 

and describes a single time step of the dynamical system. In this way a paradigmatic 
example of a dynamical system with an algebraic decay of correlation - the triangle map 
[1] was examined in reference [5]. Beside the presented method to construct a Markov 
matrix and other methods producing matrices with specific properties, there is often a 
need to construct a Markov matrix ad-hoc i.e. without incorporating any information 
about the system except the number of states d. If the construction procedure is 
a stochastic process then the resulting matrix is called the random Markov matrix 
and the set of such matrices form the ensemble of random Markov matrices. These 
matrices are usually used, without much theoretical background, for testing purposes. 
For example testing of algorithms or certain statistical hypothesises, where application 
to the field of dynamical systems, connected to ergodicity and mixing property, are the 
most interesting to us. In the information theory random Markov matrices are used to 
test the algorithms of recognition or attribution process, compression algorithms etc. 

The presented work is strongly related to the work of [6] discussing the general 
properties of the Markov ensemble and to references [7j and |8j , where a closer look at 
the second largest (sub-dominant) eigenvalue of random Markov matrices was taken. 
In contrast to the past works we try to emphasise the physical application of results in 
particular in the field of dynamical systems. Interestingly, recently ensemble of random 
Markov matrices was applied in the spectral analysis of the random quantum super- 
operators [9]. 

2. Preliminaries 

The set of all possible Markov matrices M of dimension d is defined as 

M{d) = {M e Rf"^ : Ml = 1} , 1 = (1, . . . , 1) G 
and it is isomorphic to the direct product of d convex sets 

d 

G M+ : x^l = 1} . 

The set Ai{d) forms together with the matrix product a semi-group, whereas the set 
of non-singular Markov matrices form a group of stochastic matrices. The ensemble 
of random Markov matrices is defined as a set of Markov matrices Ai{d) with the 
probability measure of matrix elements Mj j G M+ reading 

d 

dP{M) = [{d - 1)!]'^5^(M1 - l)dM , dM := JJ dMij , (1) 
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which incorporates minimal information about the set i.e. only constrains due to the 
probability conservation in the Markov process. The ensemble of random Markov 
matrices, denoted by the pair {dP{M), Ai{d)), is also referred to as the Dirichlet 
ensemble [6] and corresponding matrices are called doubly (row and column) stochastic 
or bi-stochastic Markov matrices [2\. 

The rows of the Markov matrix from the ensemble {dP{M), Ai{d)) are independent 
random vectors X = (Xj > 0)^^^^ with the distribution 

Prows(X) = (rf-l)!5(l^X-l). (2) 

It can be rather awkward to numerically generate components of vector-rows X directly 
and so a different approach to do that is taken. By following the notes on the exponential 
distribution in [in], p. 76, we find that the vector- rows X of Markov matrices can be 
expressed by vectors Y = (i^j)f=i of independent variables with a common exponential 
distribution in the following way 

^ = 1^- (3) 

This the way we numerically generate the pseudo-random rows in a Markov matrices, 
where each row generated independently. Consequently, the distribution of the rows can 
be written as 

Prows(X) = / d'YS" (x - e-i"^ . (4) 

This identity of expressions is checked by calculating the moment-generating function 
the expression above and expression ^ yielding the same result equal to 

/ d'^Xe-^'^P.owslX) = {-ir-\d - 1)! ^ . 

where A = {\i)f^i and w{x) = Y['i=ii^ ~ -^j)- Let us denote the sum of variables Yi 
by 5" = l^Y and examine its statistics. The ratio between the standard deviation 
as = a/((5' — (5*))^) of the S and its average value {S) is equal to as/ {S) = d~^ and 
it is decreasing with increasing dimension d. This means that the renormalization of 
variables Y^ in the expression ([s]) has less and less effect on the functional dependence 
between X and Y, as d is increased. We conclude that in the limit for large dimensions 
d ^ 1 variables Xj are approximately independent and exponentially distributed with 
distribution Pin(Xj) = d exp(— c/Xj). Following this idea we write the asymptotic 
approximation of the probability measure of Markov ensemble as 

dP(M) ~ dPasym(M) = d^\-^^^^^^ dM , (5) 

which is valid in the limit d oo. It can be verified that the averages w.r.t. the 
distributions dP(M) and dPasym(^) of an well behaved observable on A4, which 
depends only on a finite number of matrix elements, are asymptotically equal. 
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The probability measure of the Markov ensemble dP(M) has several unique 
properties that makes the defined ensemble of random Markov matrices {dP{M),Ai{d) 
interesting and potentially useful. 

For instance the probability measure dP{M) has a maximal Shannon entropy and 
in this information sense it is unique. The set of Markov matrices is merely a direct 
product of planes restricted to M"^ and dP{M) is uniform on them. The Shannon entropy 
of the measure dP{M) on the set Ai{d) is just the sum of Shannon entropies of the 
uniform distribution on the planes, which are themselves minimal. Hence the Shannon 
entropy of dP{M) is also minimal. Any modification of the measure would necessarily 
increase the Shannon entropy and therefore it is unique. 

It is also interesting to notice that the measure dP{M) is not invariant w.r.t. matrix 
multiplication. However, for a given non-singular Markov matrix A ^ Ai{d) the measure 
dP{M) is invariant on matrix multiplication up a constant 

P{AB) = I det{A)\-'^P{B) C M{d) . 

In fact there is no measure of Markov matrices with the matrix elements Mjj 
approximately independent in the limit d ^ 1, which would be invariant w.r.t. matrix 
multiplication. To show this let us consider two large matrices A = [^ij]fj=i and 
B = with the matrix elements being i.i.d. variables with the distribution 

P^{x). Here we denote i-th central moment of some distribution Q{x) for i > 1 as 
f^i{Q) = J dx {x — fii{Q)yQ{x) and for i = 1 by fii{Q) = J dxQ{x)x. We assume that 
the first three central moments of Pm{x) are finite. The matrix elements of the product 
AB = [Cij = X]fc=i '^i,kBk,j]fj=i are distributed by 

f (i(2 / x\ 

Pc{x) = Pab * ■ ■ ■ * Pab{x) , Pj^^{x) = — Pm(a)Pm - da, 
^ V ' a \aJ 

d + 

where the sign * denotes the convolution and Pab is the distribution of the product 
AijCj^k with the first two central moments reading: 

/^i(Pab) = f^liPraf and fl2{PAB) = 2 (/il(Pm)' + /i2(Pm)) /i2(Pm) • 

According to the central limit theorem (CLT) [10] the distribution Pc{x) converges in 
the limit oo to a Gaussian distribution with the first two central moments equal 
to yUi(Pc) = (i/ii(PAB) and /i2(Pc) = dfi2{PAB)- For distribution Pm{x) to be invariant 
w.r.t. to the matrix multiplication it has to be asymptotically equal to Pc{x) meaning 
that Pm{x) is also a Gaussian distribution for large dimension d ^ 1. By comparing the 
first two central moments of Pc{x) and Pm{x) we find that the average value of matrix 
elements of the Markov matrix and their variance are asymptotically equivalent to 

1 12 

/il(Pm) ~ 2 Ai2(Pm) 2~ ~d? ' 

respectively. The ratio between the standard deviation and the average scales with 
dimension as \/ ^i2{Pm)/ l^i{Pm) = 0{d^) and diverges in the limit d ^ oo. This indicates 
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that a measure of Markov matrices by which the matrix elements are asymptotically 
independent and distributed by Pm{x) does not exist. 

In the following we discuss the properties of random Markov matrices from the 
Markov ensemble {dP{M), A4). We focus on the entropy growth rate and correlation 
decay in the Markov chains generated by these Markov matrices, and examine their 
asymptotic behaviour for c? 3> 1. 

3. The entropy and stationary distribution of the random Markov matrices 

We consider a Markov chain defined on the set of states S = {si}f^i and with the 
conditional probabilities P{sj\si) = Mij given in the Markov matrix M = [Mi j]fj^^. 
The initial probability distribution over the states is {P{si))f^i. The probability that 
the Markov chain has evolved up to time t following a specific route (ei, . . . , Cj) G 5* is 
given with the product of conditional probabilities reading 

P(ei, ...,et)= P(ei)P(e2|ei)P(e3|e2) . • . P(ei|et_i) . 

Then the dynamic entropy S of the Markov chain at time t is given by the sum 

S{t) = -J2Pie) logP(e), 

taken over all different routes up to time t. In ergodic Markov chains we expect that 
the entropy in the limit t —>■ oo increases linearly with increasing time t as 

S{t) ~ ht, 

where G M denotes the asymptotic entropy growth rate of the Markov chain. The 
entropy growth rate is given by the formula [TT] 

d d 

/i = -^7ri5^M,,,logM,,,, (6) 
i=i j=i 

where we use the stationary distribution vr = (tTj > 0)^^^ defined as the eigenvector of 
the Markov matrix corresponding to the unit eigenvalue, 

d 

vr^M = 71^ , ^ = 1 . (7) 

1=1 

In the following we discuss the distribution P.jr{x) of elements tTj corresponding to a 
stationary distribution vr of a random Markov matrix M = [Mij]'lj^-^. In particular 
we are interested in the asymptotic limit as ^ oo, where the matrix elements Mij 
are approximately independent variables with an exponential distribution P^{Mi j) = 
dexp(—dMij). Further we assume that pi and Mij have no correlations. By doing this 
the eigenvalue equation (j?]) written by components vTj = vr^Mj j can be translated 
into an asymptotic self-consistent condition for the distribution Pn{x) reading 

P^(x) ^ P^M * . . . * Pnu ix) as d-^oo, (8) 

^ V ' 

d 
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with the distribution PttmI^;) 
PTrix) and written as 



of the products PjMj^i depending again on distribution 



Pnuix) = d J y exp y-dj^J Pn{b) , 

where * denotes the convolution. Assuming that the first three central moments of the 
distribution Pt,m{x) are finite we can use the CLT and state that the distribution Pn{x) 
converges to a Gaussian distribution as c? — oo. By inserting the ansatz 



PM 



exp 



[X — TT) 

2(yl 



n 



1 



into equation ([8j) and imposing the asymptotic equality ([8j) we obtain the variance 
o"^ ~ d^^ of the elements vTj. By appropriately rescaling the coefficients tTj their 
cumulative distribution is independent of dimension d in the limit d ^ oo and reads 



Prob 



dTCi — 1 



< 



G{x) 



X 



+ 1 



(9) 



V2 Vd 

This result is compared in figure [T] with numerically obtained distributions of rescaled 
coefficients tTj for different large dimensions d and we find a very good agreement. We 
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Figure 1. The cumulative distribution of the rescaled coefficients tt^ corresponding to 
a stationary distributions tt — {Tri)f^i (a) and its deviation form the expected limiting 
form G{x) ^ (b) calculated for individual random Markov matrices of different 
dimensions d taken from the ensemble. 

continue the analysis of the entropy growth rate /i (|6]) of an typical random Markov 
matrix from the ensemble by decomposing it into an average term /lave and an oscillating 
term hose reading 

^ d d . 

h /iave ~l~ ^osc i ^ave ^ ^ ^ Ui , /^osc ^ ^ f ^ ) iUi ^ave) ; 

i=l i=l ^ ^ 

where we introduce auxiliary variables 

d 

?7i = -^MijlogM,,,-. 
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In the asymptotic limit d ^ oo the variables Ui have according to the CLT a Gaussian 
distribution with the first two central moments reading 

m ~ log(e^-M) , 

al, ~ [1 + (7 - 4)7 + ttVS + (27 - 4 + log d) log d]^ = 0{d-^) , 

with 7 being the Euler constant. The average (■) in the expressions above is taken 
w.r.t. the asymptotic distribution Pasym(^) It is easy to see that the average term 
converges with increasing d to (Ui) as 

/iave= (f/i)+0(rf-5), 

where the last term on the r.h.s. denotes the statistical deviation. The oscillating term 
/lose can be treated as a scalar product of vectors (tTj — l/d)f^i and {Ui — /iosc)f=i and 
by applying the Schwarz-Cauchy inequality it be bounded from above: 

The last term on r.h.s. denotes again the statistical deviation obtained by taking into 
account the statistical deviations {{ni — l/d)"^) = 0{d~~^) and {{Ui — /lave)^) = 0{d~^). 
We conclude from this analysis that h behaves over the ensemble, in the leading order 
in d, approximately as the Gaussian variable /lave or Ui. This is supported by the 
numerical results in figure [2ja, where we show that the cumulative distribution of 
{h — {h))/ah converges to the cumulative normal distribution G{x) By averaging h 
over the ensemble the fluctuations from {Ui) are eliminated and we obtain the asymptotic 
equivalence 

{h) ~ log{e'^^^d) as d ^ 00 . 

The later agrees very well with the numerically obtained {h) as we can see in figure 
[2jb. This result supports the rule of thumb saying that in order to describe a dynamical 
system with an entropy rate h accurately via a Markov process we need to a Markov 
matrix of the size d ^ e^. The standard deviations cr^ scales as 0{d~^), as seen in figure 
[2]c, and this agrees with the expected statistical deviation of the form \/(ofrT~/^l^j7^ 
in our asymptotic approximation. 

4. The correlation decay induced by random Markov matrices 

A state of a Markov chain defined by a Markov matrix M = [Mij]fj^-^ G is 
described by a probability distribution 

d 
i=l 




Figure 2. The cumulative distribution of the entropy growth rate h of random Markov 
matrices for different dimensions d (a), and the average entropy growth rate (h) (b) 
and its standard deviation an (c) as a function of dimension d calculated using TV — 10^ 
random Markov matrices from the ensemble. 



over given set of states {si}f^^. Some initial probability distribution p G M^^ is evolved 
in time to p{t) by the Markov matrix in the following way 

p{tf = p^M' , 

where t G Nq denotes the discrete time. We find that a Markov chain generated by a 
typical random Markov matrix M is mixing and consequently ergodic [1]. We assume 
that the measure of Markov matrices in the ensemble corresponding to a non-mixing 
Markov chain is zero. 

The discrete analogue of the time correlation function Cf^g{t) between two real 
observables f = {fi E ^)f=i and g = {Qi & I^)iLi is defined as 

CfAt) = {gim, - {g^), , 

where we introduce in time propagated observable git)'^ = g'^M^ and averaging over 
the stationary distribution (wj). = Ylii^^i'^i- The second largest eigenvalue (called also 
the sub-dominant eigenvalue) of the Markov matrix z/ G C determines the decay of 
correlation between almost all pairs of observables (/, g) following the formula 

|Cy,g(t)| =0(|z/|*) =0(e-*/"^) as t ^ oo , 
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with Tc = — log|z/| called the correlation decay time. It is important to notice that 
spectrum A = {A : det(M — Aid.) = 0} of a Markov matrix M has the symmetry 

A* = A, 

where (■)* represents the complex conjugation. The symmetry can be noticed in figure 
(|3| where we show a spectrum of a typical random Markov matrix in the complex plane. 
In the limit of large dimensions c? 3> 1 the eigenvalues are distributed symmetrically 
around the origin with the constant distribution of the square absolute value of the form 

Prob(x < ||A|p < a; + da;) ^ 0{d~'^)dx . 

This feature goes along the lines of the Girko's circular law [12] and its generalizations 
[13] . but this particular case is not yet proved to the best of our knowledge. Here we 
are mainly interested in the second largest eigenvalues G C of the random Markov 
matrices. These are depicted for = 10^ matrices sampled uniformly across the 
ensemble in figure |3jb for several different dimensions d. For large d the values of 
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Figure 3. The spectrum of a random Markov matrix of dimension d = 10^ (a) without 

the eigenvalue 1, which corresponds to the stationary distribution, and the second 
largest eigenvalues ly E C (b) of approximately TV = 10^ random Markov matrices 
calculated for different dimensions d. 



V are distributed radially symmetrically around the origin with the average radius 
and dispersion decreasing with increasing d. Further we examine the distribution of 
the magnitude of the second largest eigenvalue \v\ denoted by P\y\ and its first two 
central moments: average magnitude (|z/|) and standard deviation o\y\. The cumulative 
distribution of the rescaled magnitude 

, _ IH - M 

0\v\ 

is depicted in figure |4}a. From the figure we conclude that the distribution of ^ is 
basically independent of dimension d for large d and we find that it agrees well with the 
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extreme value statistics of type 1 (Gumbel) [H]. Let us assume Xi are i.i.d. standard 
Gaussian variables. Then the maximal value of d variables y = max{a;j}f^^ is distributed 
according to the cumulative distribution 

where we use cumulative Gaussian distribution G{y) (|9|. It is known that under 
simple linear transformation of the variable y, which depends on d, the distribution 
of transformed y converges in the limit d ^ oo to the Gumbel or double exponential 
distribution. To avoid certain slow convergence problems towards the limiting 
distribution outlined in [15J we compare in figure |4ja numerically obtained distribution 
directly with Pmax((y — y)/<^y,d) for several large enough d, where y and ay are the 
average maximal value and its standard deviation, respectively. We find a very good 
agreement suggesting that the eigevalues of the Markov matrix behave as i.i.d random 
complex variables inside some disk in the complex plane with the radius 0{d~2). The 
first two central moments of numerical results as function of d are shown in figures |4jb 
andlllc. The average magnitude of the second largest eigenvalue fits very well to 
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Figure 4. The cumulative distribution of the magnitude of the second largest 
eigenvalue \i>\ of random Markov matrices at different dimensions d (a) and the average 
(Ivl) (b) and the standard deviation of the amplitude (T|y| (b) as a function of the 
dimension d calculated from approximately 10^ and 10"^ random Markov matrices from 
the ensemble for d < 512 and d > 512, respectively. 



the empirically found asymptotic formula 



(|z/|)~Corf- 



d 



oo 



(10) 
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where Co ~ 1. This resuh supports the conjecture that the Girko's circular law is valid 
for the random Markov matrices stated in |6j and it can be better understood by the 
asymptotics of the upper bound for |z/| obtaining in following. By taking into account 
that all left-hand eigenvector except the stationary distribution n are perpendicular to 
the vector 1, we can upper bound of the second largest eigenvalue as 

< limsup \\x'^M\\l , S = {x eC^ : \\x\\2 = 1 A x ± 1} . 

Here we write the expression ||a;"'^M||| = x^Nx using the matrix N = MM'^. In the 
asymptotic limit the matrix takes the form 

d 

k=l 

where the last term denotes the statistical error of the expression. From here we 
immediately obtain asymptotic expression of the upper bound 

limsup ||Mx||2 ~ . 

This means that the second largest eigenvalue in a typical random Markov matrix is 
bounded by below d~^ in the limit d oo. This is true also in the average over the 
ensemble yielding the relation (z/) < d^^ and setting the value of the constant Co = 1. 
The asymptotic behaviour of the standard deviation (T\^\ is not as clear as in the case of 
the average value The numerical results suggest the power law decay 

cx\u\ ~ Ci , a ^ 1.1 . 

In the Markov approximations of a dynamical systems with the mixing property [TE] , as 
described in the introduction, it interesting to know about the correlation decay-time Tc 
and entropy growth rate h and their dependence on the cardinality of the state space. 
In the random Markov matrices form the ensemble we find that the average average 
correlation decay-time t := — l/log(|z/|) and average entropy growth rate {h) obey the 
asymptotics 

(/;,) Tc ~ - as d ^ oo . 

In dynamical systems there are strong indications that the correlation decay time and the 
entropy growth rate, given as the sum of positive Lyapunov exponents, are correlated, 
but this connection is not well understood, yet, see e.g. [HI [18]. We address this 
question in the random Markov matrices and calculate pairs corresponding to 

the Markov matrices sampled uniformly over the ensemble. The result is graphically 
depicted in figure [5} where one can notice that in particular at small dimensions d there 
is clearly some correlation between the amplitude of the second largest eigenvalue 
and entropy growth rate h of random Markov matrices. 
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Figure 5. The amplitude of the second largest eigenvalue and corresponding 
entropy h calculated for some number of random Markov matrices samples uniformly 
with respect to the measure dP{M) (ij) at different dimensions d. For details on 
statistics see the caption under figure [4 



The latter is tested by calculating the statistical correlation between the reciprocal 
correlation decay time t^^ = — log |z/| and the entropy growth rate h over the ensemble 
of random Markov matrices and is given by 



Corr((i) 



{r.-')f){{h-{h)f) 



The correlation Corr((i) as a function of dimension d is presented in figure [6] and we see 
that it slowly decreases with increasing d. Currently it is not possible to determine the 
dependence of the correlation Corr((i) on dimension more precisely. 
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Figure 6. The normalized correlation Corr((i)/Corr(0) between the reciprocal 
correlation time t^^ = — log|z^| and the entropy growth rate h in the ensemble of 
random Markov matrices as a function dimension d. For details on statistics see the 
caption under figure |4] 



5. Conclusions 



We define the ensemble of random Markov matrices, present its basic properties and 
point out few of the potential physical, technical and mathematical applications. We 
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analyse the statistical properties of the stationary distribution tt = (vrj)f^^ corresponding 
to a typical element of the ensemble, and study the distribution of the entropy growth 
rate h over the ensemble, where we find a good agreement with analytical predictions 
stating that tTj is a Gaussian variable and h is asymptotically equal to \og{e'^~^d) in the 
limit of large dimensions d ^ oo. Further we analyse the second largest eigenvalue u of 
the Markov matrices, which is connected to the correlation decay in the Markov chains. 
We show numerically and provide a heuristic proof that in average over the ensemble 
the second largest eigenvalue decreases with increasing dimension d as ~ c?^5 . 
Additionally we calculate the correlation between the correlation decay rate and the 
entropy growth rate and find that it decreases with increasing dimension of the Markov 
matrices. 

We believe that the current results enrich the understanding of Markov processes in 
the limit of large state spaces and all applications which can be described by Markov 
processes. 
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